Ejercicio 1

Discuta con sus compañeros, ¿Qué fenómenos físicos pueden influir en la relación que estamos estudiando entre estas dos variables? ¿Son lógicos los resultados de estas mediciones?

Motivación de la regresión lineal

Uno de los objetivos de la ciencia consiste en predecir y describir sucesos del mundo en que vivimos. Una manera de hacerlo es construir modelos matemáticos que describan adecuadamente el mundo real, y la regresión lineal no es la excepción. Si analizamos con detenimiento la gráfica de dispersión de los datos de Presion y punto de Rocio podemos establecer con claridad una relación entre estas dos variables. Si \(y\) representa el punto de rocio y \(x\) la presion para algún municipio en la región Orinoquía de Colombia, podemos crear un modelo determinístico de la siguiente forma:

\[y = \beta_0 + \beta_1x\] De hecho, da la impresión que los datos caen, en general, pero no exactamente, en una línea recta, por lo cuál, deberíamos modificar dicha ecuación para tener en cuenta el error aleatorio que rige todo el estudio de estos datos (puede consultar todos los tipos de errores en Holmes, Tyson (2004)). Por lo tanto un modelo que se aproxima más a la realidad sería un modelo probabilístico

\[Y = \beta_0 + \beta_1x + \epsilon\]

La ecuación anterior se llama Modelo de Regresión Lineal Simple. Analicemos con detenimiento cada parte de nuestra ecuación:

• \(Y\) variable dependiente
• \(\beta_0\) y \(\beta_1\) coeficientes de regresion: Valores fijos
• \(x\) variable independiente: Valor fijo conocido, podemos medirlo.

La componente de error que acompaña la ecuación es aleatoria y determina las propiedades de nuestra variable respuesta. Supongamos que \(\epsilon \sim (0, \sigma^2)\). Por lo tanto la respuesta esperada de \(y\) dado \(x\) es:

\[E(y|x) = \mu_{y|x} = E(\beta_0 + \beta_1x+\epsilon) = \beta_0 + \beta_1x\] Así, el verdadero modelo de regresión \(\mu_{y|x} = \beta_0 + \beta_1x\) es una línea recta de valores promedios, esto es, la altura de la línea de regresión en cualquier valor de \(x\) no es más que el valor esperado de \(y\) para esa \(x\). Se puede interpretar que la pendiente \(\beta_1\) es el cambio de la media de \(y\) para un cambio unitario de \(x\). Además, la variabilidad de \(y\) en algún valor particular de \(x\) queda determinada por la varianza del componente de error en el modelo, \(\sigma^2\). Esto implica que hay una distribución de valores de \(y\) en cada \(x\), y que la varianza de esta distribución es igual en cada x.

Cuando \(\sigma^2\) es pequeña, los valores observados del tiempo de entrega serán cercanos a la recta, y cuando \(\sigma^2\) es grande, se pueden desviar bastante de la línea.

Explicación modelo de regresión lineal simple

Para que el modelo de regresión lineal simple sea válido, debe cumplir 4 supuestos:

1. Relación lineal: La relación entre la variable dependiente y la independiente debe ser lineal. Esto puede apreciarse con algo tan sencillo como un diagrama de puntos o de dispersión, que nos muestra el aspecto de la relación en el rango de valores observados de la variable independiente. Sin embargo, cada supuesto debe examinarse a través de métodos estadísticos de manera formal.

2. Normalidad del error: La distribución aproximada de los errores, debe ser normal.

3. Errores no autocorrelacionados: Los residuos deben ser independientes entre sí y que no haya ningún tipo de correlación entre ellos.

4. Homocedasticidad del los errores: Los residuos deben distribuirse de forma homogénea para todos los valores de la variable de predicción.

Después de obtener el ajuste por mínimos cuadrados, surgen varias preguntas interesantes:

1. ¿Qué tan bien se ajusta esta ecuación a los datos?

2. ¿Es probable que el modelo sea útil como predictor?

3. ¿Se viola alguna de las hipótesis básicas (como la de varianza constante y (como la de varianza constante y la de errores no correlacionados)? y en caso afirmativo, ¿qué tan grave es eso?

Los residuales juegan un papel clave para evaluar la adecuación del modelo. Se puede considerar que los residuales son realizaciones de los errores \(\epsilon_i\) del modelo. Así, para comprobar la constancia de la varianza y la hipótesis de errores no correlacionados, uno se debe preguntar si los residuales parecen ser realmente una muestra aleatoria de una distribución con esas propiedades.

Identificación

Volvamos a los datos de Presion y Punto de Rocío para diferentes municipios en la región Orinoquía. Para la identificación del mismo, es necesario realizar un análisis descriptivo completo, por esta vez omitimos este análisis, ya que previamente se ha realizado en el Taller I y al inicio de este notebook.

library(ggrepel)
ggplot(data4, aes(x=Presion, y=Punto_de_Rocio)) +
  geom_point() + 
  geom_text_repel(aes(Presion, Punto_de_Rocio, label = rownames(data4)),max.overlaps = 12,
                      size         = 2,
                      box.padding  = 0.5) 

## Warning: ggrepel: 14 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

Como vimos anteriormente, existe una relación lineal entre las dos variables de análisis, por lo tanto podríamos aplicar un modelo de regresión lineal simple.

Estimación

Estimación desde R

Para la estimación del modelo ocuparemos la función lm() que se encuentra incorporada entre las funciones base de R.

ggplot(data4, aes(x=Presion, y=Punto_de_Rocio)) + 
  geom_point()+
  geom_smooth(method=lm, se=FALSE)

## `geom_smooth()` using formula = 'y ~ x'

reg <- lm(Punto_de_Rocio ~ Presion, data = data4)
summary(reg)

## 
## Call:
## lm(formula = Punto_de_Rocio ~ Presion, data = data4)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1279 -0.4386 -0.0305  0.5099  1.2279 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1319.70436   63.31502   20.84   <2e-16 ***
## Presion       -1.28251    0.06255  -20.50   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6308 on 57 degrees of freedom
## Multiple R-squared:  0.8806, Adjusted R-squared:  0.8785 
## F-statistic: 420.5 on 1 and 57 DF,  p-value: < 2.2e-16

Salida de la estimación del modelo en R

1. Identificación del modelo

2. Estadísticas descriptivas de los errores calculados con el modelo

3. Tabla de coeficientes asociados al modelo de regresión lineal, por las filas encontramos todos los parámetros asociados al modelo. Por las columnas encontramos:

• Estimación: Valor de la estimación asociado al parámetro correspondiente.

• Std. Error: Error estándar de la estimación, se usa para construir intervalos de confianza para este parámetro.

• t-value: Estadístico de prueba t, utilizado para evaluar la hipótesis:

\[H_0 : \beta_i = 0\] \[vs\] \[H_a : \beta_i \neq 0\]

La forma de interpretar los coeficientes es la siguiente:

1. Intercepto: Con un p-valor superior a 0.05, no podemos rechazar la hipótesis nula. Por lo tanto se asumir que este no es importante en el modelo.

2. Presion: Con un p-valor inferior a 0.05, concluimos que existe suficiente evidencia estadística para rechazar \(H_0\). Con lo cuál, interpretamos que la Presion explica el punto de rocio y es importante dentro de nuestro modelo.

• Pr(>|t|): P-valor, asociado al estadístico de prueba incluido en el sistema de hipótesis anterior.

4. Resultados para el ANOVA asociado al modelo. El análisis de varianza permite determinar la bondad de ajuste del modelo a los datos.

anova(reg)

## Analysis of Variance Table
## 
## Response: Punto_de_Rocio
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## Presion    1 167.308 167.308  420.46 < 2.2e-16 ***
## Residuals 57  22.681   0.398                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Del ANOVA asociado al modelo podemos extraer varias conclusiones:

• El error estándar residual viene dado por \(\hat{\sigma} = \sqrt{\dfrac{\sum \hat{\epsilon}^2}{n-p}}= \sqrt{\dfrac{\sum y_i -\hat{y}_i^2}{n-p}}\). Da la desviación estándar de los residuos y nos dice qué tan grande es el error de predicción en la muestra o en los datos de entrenamiento.

• \(R^2\): Descrito a detalle anteriormente.

• \(R^2\) ajustado: Note que incorporar una variable al modelo, esta puede o no ser relevante para explicar la variabilidad en \(y\). Si la variable \(z\) no tiene efecto en la respuesta:

• Al minimizar \(SC_M\), se lleva a \(\beta_2=0\) y así \(y = \beta_0 + \beta_1 x+ 0 z \rightarrow y = \beta_0 + \beta_1 x\) tenemos el modelo de RLS.

• \(SC_M\) es la misma bajo los dos modelos, es decir añadir \(Z\) no mejoró, ni empeoró \(R^2\).

• Añadir variables mantiene igual o incluso mejora \(R^2\) aunque no sean de utilidad para explicar \(Y\).

  • Variables irrelevantes pueden estar correlacionadas con la variable \(Y\) por coincidencia.
  • Mas variables irrelevantes aumentan la probabilidad de que esto suceda, mejorando artificialmente \(R^2\).

En la práctica se reporta usualmente un \(R^2\) ajustado por en número de variables que está dado por:

\[R^2_{adj} = 1- (1-R^2)\dfrac{n-1}{n-p-1}\]

La interpretación es la misma que la del \(R^2\)

• El estadístico de prueba F asociado al modelo es el encargado de evaluar la hipótesis:

\[H_0 : \beta_0 = \beta_1 =...=\beta_p = 0\] \[vs\] \[H_a : Existe\ al\ menos\ un\ i\ tal\ que\ \beta_i \neq 0\] Lo cuál evalúa en conjunto todo el modelo, en este caso en particular como el P-valor asociado a la prueba es menor que un nivel de significancia de 0.05, podemos rechazar la hipótesis nula de la prueba y por lo tanto asumir que el modelo está bien ajustado a los datos.

Validación

Para la validación del modelo, debemos aplicar diferentes estrategias según sea el caso del supuesto que estemos evaluando:

par(mfrow=c(2,2))
plot(reg)

par(mfrow=c(1,1))

1. Del gráfico de Residuales vs valores ajustados, esperamos patrones aleatorios, esto es sinónimo de homocedasticidad o varianza constante en los errores.

2. Del gráfico Q-Q plot esperamos un comportamiento uno a uno sobre la recta, este gráfico compara percentiles teóricos de una distribución normal, con los obtenidos en los errores por el modelo, si logramos obtener que todos los puntos están sobre la recta, podemos asumir que cumplimos el supuesto de normalidad del error.

3. Del gráfico Scale-Location obtenemos una conclusión similar al primero, que no haya tendencia en el gráfico (una línea horizontal sin pendiente) y que los puntos alrededor de esta línea no presenten patrón alguno.

4. Del último gráfico podemos hacer un análisis de influencias, puntos que se encuentren en la parte derecha de este gráfico, son probablemente datos que podrían estar afectando la calidad del modelo, y debería contemplarse la posibilidad de excluir estos puntos del modelo. Esto se logra, comparando los modelos con y sin estos posibles datos influyentes.

5. Finalmente para evaluar la autocorrelación de los errores de forma gráfica, el gráfico de errores vs ajustados de la observaciones debe mantener el mismo componente aleatorio que hemos visto en los supuestos anteriores.

Uso

Finalmente, cuando el modelo ha sido probado y ha cumplido con todos los supuestos, solo en ese momento podemos hacer uso del mismo.

Suponga que se tiene una predicción de la presión para alguno de los municipios para el día siguiente con un valor de 1013.5; cuál sería el valor esperado de Punto de Rocío dado este valor de presión.

predict(reg, newdata = data.frame(Presion = 1013.5))

##        1 
## 19.87916

Por lo tanto el punto de rocío dado ese valor en la presión para algún municipio de la región Orinoquia en particular es de 19.8

Ejercicio 2

Ejecuta una regresión lineal múltiple tratando de explicar el punto de rocío en función de la presión y la temperatura.

Como se puede observar, la mejora en el ajuste del modelo al agregar una nueva variable es bastante claro, mejoramos tanto el \(R^2\) como el \(R^2_{adj}\); con lo cuál podemos pensar que el modelo logró ajustrarse mejor a la nuestra variable respuesta.

Para entender un poco mejor como funciona en un espacio de 3 dimensiones este nuevo modelo de regresión lineal, graficaremos a continuación las estimaciones como un plano proyectado en este nuevo espacio:

mesh_size <- .5
margin    <- 0

x_min <- min(data5$Presion) - margin
x_max <- max(data5$Presion) - margin
y_min <- min(data5$Temperatura) - margin
y_max <- max(data5$Temperatura) - margin
xrange <- seq(x_min, x_max, mesh_size)
yrange <- seq(y_min, y_max, mesh_size)
xy <- meshgrid(x = xrange, y = yrange)
xx <- xy$X
yy <- xy$Y
dim_val <- dim(xx)
xx1 <- matrix(xx, length(xx), 1)
yy1 <- matrix(yy, length(yy), 1)
final <- data.frame(cbind(xx1, yy1))
colnames(final)<- c("Presion", "Temperatura")

pred <- mod %>%
  predict(final)

pred <- matrix(pred, dim_val[1], dim_val[2])

plot_ly(data5, x = ~Presion, y = ~Temperatura, z = ~Punto_de_Rocio ) %>% 
  add_markers(size = 5) %>% 
  add_surface(x=xrange, y=yrange, z=pred, alpha = 0.65, type = 'mesh3d', name = 'pred_surface')

Ejercicio 3

Verifique los nuevos supuestos de este modelo de regresión con nuevas variables explicativas.

Finalmente para evaluar la multicolinealidad en el modelo podemos acudir a dos metodologías que podrían dar indicios de estos problemas:

1. Matriz de Correlación: Encontrar correlaciones lineales con valores cercanos a un valor de 1 (o -1) puede dar indicios de problemas de multicolinealidad en nuestro modelo:

cor(data5 %>% dplyr::select(-Punto_de_Rocio))

##                Presion Temperatura
## Presion      1.0000000  -0.9552992
## Temperatura -0.9552992   1.0000000

2. Factores de inflación de varianza: Esta metodología permite detectar multicolinealidad en nuestro modelo determinando si las varianzas de los coeficientes de regresión están infladas por problemas de multicolinealidad. Los detalles técnicos de esta metodología se pueden consultar en Montgomery, D., Peck, E. A., & Vining, G. G. (2006).

vif(mod)

##     Presion Temperatura 
##    11.44119    11.44119

Existe una regla empírica que establece que valores superiores a 10 en este indicador, demuestran problemas de multicolinealidad. El paso a seguir, debería ser extraer dichas variables o correr un modelo de regresión tipo Ridge o Lasso que puedan resolver de otra forma este conflicto entre variables explicativas.